;+
;  Creates a Mandelbrot fractal object that can be iterated, zoomed and
;  displayed in either a direct graphics or object window. See
;  `this MathWorld page <http://mathworld.wolfram.com/MandelbrotSet.html>`
;  for further details on the Mandelbrot set.
;
; :Categories:
;    Graphics
;
; :Examples:
;   Direct graphics window version::
;     imageSize = 300
;
;     theFractal = OBJ_NEW('DJP_Mandelbrot', imageSize=imageSize)
;
;     ;; Iterate the fractal 5 times then draw it.
;     theFractal -> Iterate, 5
;     cgLoadCT, 5, /BREWER, /REVERSE
;     WINDOW, 0, XSIZE=imageSize, YSIZE=imageSize
;     theFractal -> Draw
;
;     ;; Iterate the fractal another 25 times (30 altogether) and
;     ;; then draw it in a new window
;     theFractal -> Iterate, 25
;     WINDOW, 1, XSIZE=imageSize, YSIZE=imageSize
;     theFractal -> Draw
;
;     OBJ_DESTROY, theFractal
;
;   Object graphics version::
;     imageSize = [300, 300]
;
;     oDraw    = OBJ_NEW('IDLgrWindow', DIMENSIONS=imageSize)
;     oDrawRef = PTR_NEW(oDraw)
;
;     oView    = OBJ_NEW('IDLgrView', VIEWPLANE_RECT=[0, 0, imageSize[0], imageSize[1]])
;     oViewRef = PTR_NEW(oView)
;
;     oModel   = OBJ_NEW('IDLgrModel')
;     oPalette = OBJ_NEW('IDLgrPalette')
;
;     brewerfile = Find_Resource_File('fsc_brewer.tbl')
;     IF brewerfile EQ '' THEN oPalette -> LoadCT, 23
;     IF brewerfile NE '' THEN oPalette -> LoadCT, 23, FILE=brewerfile
;
;     oImage = OBJ_NEW('IDLgrImage', PALETTE=oPalette)
;     oImageRef = PTR_NEW(oImage)
;
;     oModel -> Add, oImage
;     oView  -> Add, oModel
;
;     theFractal = OBJ_NEW('DJP_Mandelbrot', imageSize=imageSize)
;     theFractal -> Iterate, 100
;     theFractal -> Draw, oImageRef=oImageRef, oWindowRef=oDrawRef, oViewRef=oViewRef
;
;     OBJ_DESTROY, theFractal
;     PTR_FREE, oImageRef, oDrawRef, oViewRef
;
; :Author:
;    David Platten
;
; :Copyright:
;     Copyright (c) 2013, David J Platten
;
; :History:
;     Change History::
;        Written, May 2013, DJP
;-


;+
; :Description:
;   This method cleans up the pointers used in the DJP_Mandelbrot class.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot::Cleanup

  PTR_FREE, self.fractal
  PTR_FREE, self.indices
  PTR_FREE, self.zed
  PTR_FREE, self.coords
  PTR_FREE, self.extra

END

;+
; :Description:
;   This method initialises the object. It sets up the `self` structure
;   to contain the fractal default settings, as shown by the keywords.
;
; :Keywords:
;    fractal: in, optional, type=UintArr(imageSize, imageSize), default=0
;       A two-dimensional array containing the number of iterations that have
;       been used for each location in coords so far.
;    indices: in, optional, type=UintArr(imageSize, imageSize)
;       A two-dimensional array containing the indices of fractal that are
;       still being worked on.
;    zed: in, optional, type=ComplexArr(imageSize, imageSize), default=0
;       The value of z at every complex coordinate included in the fractal.
;    coords: in, optional, type=ComplexArr(imageSize, imageSize)
;       The complex coordinates of the area being calculated.
;    imageSize: in, optional, type=integer, default=500
;       The number of elements in each row and column of the fractal array.
;    rMin: in, optional, type=double, default=-2.0
;       The minimum real value to be considered.
;    rMax: in, optional, type=double, default=0.5
;       The maximum real value to be considered.
;    iMin: in, optional, type=double, default=-1.25
;       The minimum imaginary value to be considered.
;    iMax: in, optional, type=double, default=1.25
;       The maximum imaginary value to be considered.
;    numIter: in, optional, type=integer, default=0
;       The number of times that the fractal has been iterated.
;    _extra: in, optional
;       The "extra" keywords for the command that is being run.
;
; :Returns:
;   1 if successful, 0 if unsuccessful.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
FUNCTION DJP_Mandelbrot::Init,  $ ;; The name of the method
           fractal=fractal,     $ ;; The integer-array fractal
           indices=indices,     $ ;; The long-array of fractal indices
           zed=zed,             $ ;; The complex array of iterated values
           coords=coords,       $ ;; The complex array space that is being examined
           imageSize=imageSize, $ ;; The fractal array dimensions
           rMin=rMin,           $ ;; The minimum real value
           rMax=rMax,           $ ;; The maximum real value
           iMin=iMin,           $ ;; The minimum imaginary value
           iMax=iMax,           $ ;; The maximum imaginary value
           numIter=numIter,     $ ;; The number of iterations carried out
           _Extra=extra           ;; Holds extra keywords

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  ;; Set up an error handler
  CATCH, theError
  IF theError NE 0 THEN BEGIN
    CATCH, /Cancel
    ok = ERROR_MESSAGE(!Error_State.Msg + ' Returning...', Traceback=1, /Error)
    RETURN, 0
  ENDIF

  ;; Check that the positional parameters and keywords have values
  ;; I check the dimensions of imageSize, and make a 2 dimensional variable
  ;; imageDims from imageSize. I don't simply over-write imageSize with the
  ;; 2d variable as I don't want to change it's dimensions in the routine
  ;; that called it.
  IF N_ELEMENTS(imageSize) EQ 0 THEN imageDims = [200, 200]
  IF N_ELEMENTS(imageSize) EQ 1 THEN imageDims = [imageSize[0], imageSize[0]]
  IF N_ELEMENTS(imageSize) EQ 2 THEN imageDims = imageSize
  IF N_ELEMENTS(rMin)      EQ 0 THEN rMin      = -2.0D
  IF N_ELEMENTS(rMax)      EQ 0 THEN rMax      = 0.5D
  IF N_ELEMENTS(iMin)      EQ 0 THEN iMin      = -1.25D
  IF N_ELEMENTS(iMax)      EQ 0 THEN iMax      = 1.25D
  IF N_ELEMENTS(numIter)   EQ 0 THEN numIter   = 0
  IF N_ELEMENTS(zed)       EQ 0 THEN zed       = DCOMPLEXARR(imageDims[0], imageDims[1])
  IF N_ELEMENTS(indices)   EQ 0 THEN indices   =     LINDGEN(imageDims[0], imageDims[1])
  IF N_ELEMENTS(fractal)   EQ 0 THEN fractal   =      INTARR(imageDims[0], imageDims[1])
  IF N_ELEMENTS(coords)    EQ 0 THEN BEGIN
    ;; Work out the increments in the real and imaginary directions
    rIncrement = (rMax - rMin) / imageDims[0]
    iIncrement = (iMax - iMin) / imageDims[1]

    ;; Create a 2D array contining the real position of each element
    rCoords = (DINDGEN(imageDims[0]) * rIncrement) + rMin
    rCoords = rCoords # (rCoords * 0 + 1)

    ;; Create a 2D array containing the imaginary position of each element
    iCoords = (DINDGEN(imageDims[1]) * iIncrement) + iMin
    iCoords = (iCoords * 0 + 1) # iCoords

    ;; Create a complex array from rCoords and iCoords
    coords = DCOMPLEX(rCoords, iCoords)
  ENDIF

  self.fractal   = PTR_NEW(fractal)
  self.indices   = PTR_NEW(indices)
  self.zed       = PTR_NEW(zed)
  self.coords    = PTR_NEW(coords)
  self.imageSize = imageDims
  self.rMin      = rMin
  self.rMax      = rMax
  self.iMin      = iMin
  self.iMax      = iMax
  self.numIter   = numIter
  self.extra     = PTR_NEW(extra)

  RETURN, 1
END

;+
; :Description:
;   This method resets the real and imaginary range of an existing
;   DJP_Mandelbrot object to default settings of rMin = -2.0D,
;   rMax = 0.5D, iMin = -1.25D and iMax = 1.25D. The aspect ratio
;   of imageSize is used to adjust these to ensure that the reset
;   fractal has the correct aspect ratio.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot::Reset

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  ;; Set up an error handler
  CATCH, theError
  IF theError NE 0 THEN BEGIN
    CATCH, /Cancel
    ok = ERROR_MESSAGE(!Error_State.Msg + ' Returning...', Traceback=1, /Error)
    RETURN
  ENDIF

  ;; Reset the r and i ranges
  aspectRatio = DOUBLE(self.imageSize[0]) / DOUBLE(self.imageSize[1])
  IF aspectRatio GT 1.0 THEN BEGIN
    rCentre   = (-2.0D + 0.5D) / 2.0
    rRange    = 2.5D * aspectRatio
    self.rMin = rCentre - (rRange / 2.0)
    self.rMax = self.rMin + rRange
    self.iMin = -1.25D
    self.iMax = 1.25D
  ENDIF ELSE BEGIN
    iCentre   = 0.0D
    iRange    = 2.5D / aspectRatio
    self.iMin = iCentre - (iRange / 2.0)
    self.iMax = self.iMin + iRange
    self.rMin = -2.0D
    self.rMax = 0.5D
  ENDELSE

  self -> ResetArrays

END

;+
; :Description:
;   This method resets the arrays in an existing DJP_Mandelbrot object to
;   the default settings set by `DJP_Mandelbrot::Reset`. It also sets the
;   number of iterations carried out back to zero.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot::ResetArrays

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  zed       = DCOMPLEXARR(self.imageSize[0], self.imageSize[1])
  indices   =     LINDGEN(self.imageSize[0], self.imageSize[1])
  fractal   =      INTARR(self.imageSize[0], self.imageSize[1])

  rIncrement = (self.rMax - self.rMin) / self.imageSize[0]
  iIncrement = (self.iMax - self.iMin) / self.imageSize[1]

  ;; Create a 2D array of dimensions self.imageSize[0], self.imageSize[1]
  ;; contining the real position of each element
  rCoords = (DINDGEN(self.imageSize[0]) * rIncrement) + self.rMin
  rCoords = rCoords # (DINDGEN(self.imageSize[1]) * 0 + 1)

  ;; Create a 2D array of dimensions self.imageSize[0], self.imageSize[1]
  ;; containing the imaginary position of each element
  iCoords = (DINDGEN(self.imageSize[1]) * iIncrement) + self.iMin
  iCoords = (DINDGEN(self.imageSize[0]) * 0 + 1) # iCoords

  ;; Create a complex array from rCoords and iCoords
  coords = DCOMPLEX(rCoords, iCoords)

  ;; Free up the previous pointers
  PTR_FREE, self.fractal
  PTR_FREE, self.indices
  PTR_FREE, self.zed
  PTR_FREE, self.coords

  ;; Assign the reset variables to the self structure
  self.fractal = PTR_NEW(fractal)
  self.indices = PTR_NEW(indices)
  self.zed     = PTR_NEW(zed)
  self.coords  = PTR_NEW(coords)
  self.numIter = 0
END


;+
; :Description:
;   This method recalculates the real and imaginary ranges of the
;   existing object and then calls the `DJP_Mandelbrot::ResetArrays`
;   method to setup the new arrays.
;
; :Params:
;    magnification: in, optional, type=double
;       The amount of zoom to apply to the current real and
;       imaginary ranges. A value greater than 1.0 will zoom in,
;       whilst a value less than 1.0 will zoom back out. If no
;       value is provided then the method exits without doing
;       anything.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot::Zoom, magnification

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  ;; Set up an error handler
  CATCH, theError
  IF theError NE 0 THEN BEGIN
    CATCH, /CANCEL
    ok = ERROR_MESSAGE(!Error_State.Msg + ' Returning...', TRACEBACK=1, /ERROR)
    RETURN
  ENDIF

  ;; If no magnification has been provided then return
  ;; without doing anything.
  IF N_ELEMENTS(magnification) EQ 0 THEN RETURN

  ;; Work out the new r and i max and min values
  rRangeCurrent = self.rMax - self.rMin
  iRangeCurrent = self.iMax - self.iMin

  rRangeNew = rRangeCurrent / magnification
  iRangeNew = iRangeCurrent / magnification

  rCentre = self.rMin + (rRangeCurrent / 2.0)
  iCentre = self.iMin + (iRangeCurrent / 2.0)

  self.rMin = rCentre - (rRangeNew / 2.0)
  self.rMax = rCentre + (rRangeNew / 2.0)

  self.iMin = iCentre - (iRangeNew / 2.0)
  self.iMax = iCentre + (iRangeNew / 2.0)

  self -> ResetArrays

END


;+
; :Description:
;   This method changes the real and imaginary ranges, as well
;   as the image size. It then calls `DJP_Mandelbrot::ResetArrays`
;   so that the new arrays can be calculated.
;
; :Keywords:
;    imageSize: in, optional, type=integer, default=self.imageSize
;       The number of elements in each row and column of the fractal array.
;    rMin: in, optional, type=double, default=self.rMin
;       The minimum real value to be considered.
;    rMax: in, optional, type=double, default=self.rMax
;       The maximum real value to be considered.
;    iMin: in, optional, type=double, default=self.iMin
;       The minimum imaginary value to be considered.
;    iMax: in, optional, type=double, default=self.iMax
;       The maximum imaginary value to be considered.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot::Change, rMin=rMin, rMax=rMax, iMin=iMin, iMax=iMax, imageSize=imageSize

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  ;; Set up an error handler
  CATCH, theError
  IF theError NE 0 THEN BEGIN
    CATCH, /CANCEL
    ok = ERROR_MESSAGE(!Error_State.Msg + ' Returning...', TRACEBACK=1, /ERROR)
    RETURN
  ENDIF

  ;; If all parameters are empty then return without doing anything.
  IF (N_ELEMENTS(rMin)      EQ 0) AND $
     (N_ELEMENTS(rMax)      EQ 0) AND $
     (N_ELEMENTS(iMin)      EQ 0) AND $
     (N_ELEMENTS(iMax)      EQ 0) AND $
     (N_ELEMENTS(imageSize) EQ 0) THEN RETURN

  self.rMin = rMin
  self.rMax = rMax
  self.iMin = iMin
  self.iMax = iMax
  IF N_ELEMENTS(imageSize) EQ 1 THEN self.imageSize = [imageSize[0], imageSize[0]]
  IF N_ELEMENTS(imageSize) EQ 2 THEN self.imageSize = imageSize

  self -> ResetArrays

END


;+
; :Description:
;   This method changes the image size. It then calls `DJP_Mandelbrot::ResetArrays`
;   so that the new arrays can be calculated.
;
; :Params:
;    imageSize: in, optional, type=integer, default=self.imageSize
;       The number of elements in each row and column of the fractal array.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot::Resize, imageSize

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  ;; Set up an error handler
  CATCH, theError
  IF theError NE 0 THEN BEGIN
    CATCH, /CANCEL
    ok = ERROR_MESSAGE(!Error_State.Msg + ' Returning...', TRACEBACK=1, /ERROR)
    RETURN
  ENDIF

  ;; If the parameter is empty then return without doing anything.
  IF N_ELEMENTS(imageSize) EQ 0 THEN RETURN
  IF N_ELEMENTS(imageSize) EQ 1 THEN self.imageSize = [imageSize[0], imageSize[0]]
  IF N_ELEMENTS(imageSize) EQ 2 THEN self.imageSize = imageSize

  self -> ResetArrays

END


;+
; :Description:
;   This method iterates the existing fractal by the number of
;   iterations passed to it.
;
; :Params:
;    numIterations: in, optional, type=integer, default=1
;       The number of iterations to carry out. If no value is
;       provided then one iteration is carried out.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot::Iterate, numIterations

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  ; Set up an error handler
  CATCH, theError
  IF theError NE 0 THEN BEGIN
    CATCH, /CANCEL
    ok = ERROR_MESSAGE(!Error_State.Msg + ' Returning...', TRACEBACK=1, /ERROR)
    RETURN
  ENDIF

  IF N_ELEMENTS(numIterations) EQ 0 THEN numIterations = 1

  ;; Carry out the required number of iterations
  FOR iteration=0, numIterations-1 DO BEGIN

	*self.zed = (*self.zed)^2 + *self.coords

	f = ABS(*self.zed) GT 2.0

    ;; Find out where the elements that are greater than maxSize are.
	h = HISTOGRAM(f, BINSIZE=1, MIN=0, MAX=1, REVERSE_INDICES=ri)
	count0 = h[0]
	count1 = h[1]

	IF count1 GT 0 THEN BEGIN

	  w0 = ri[ri[0]:ri[1]-1]
	  w1 = ri[ri[1]:ri[2]-1]

	  (*self.fractal)[(*self.indices)[w1]] = self.numIter + iteration

	  IF count1 LT 5 THEN BEGIN
	    BREAK
	  ENDIF

    ;; Reset indices, zed and coords to just contain those elements
    ;; that are still below the threshold value.
	  *self.indices = (*self.indices)[w0]
	  *self.zed     = (*self.zed)[w0]
	  *self.coords  = (*self.coords)[w0]
	ENDIF

  ENDFOR

  ;; Increment the number of iterations
  self.numIter = self.numIter + numIterations

END


;+
; :Description:
;   This method draws the current fractal in a window. If a
;   reference to an image object are passed then the fractal
;   is drawn using the supplied object references to the image,
;   view and window. If no image object reference is provided
;   then the fractal is drawn using the Coyote cgImage routine
;   in the current direct graphics window. If no window is
;   open then this command will create one.
;
; :Keywords:
;    oImageRef: in, optional, type=IDLgrImage
;       A reference to an IDLgrImage object that is to be used to
;       display the fractal data.
;    oViewRef: in, optional, type=IDLgrView
;       A reference to an IDLgrView object that includes the
;       image object.
;    oWindowRef: in, optional, type=IDLgrWindow
;       A reference to an IDLgrWindow object that is to display
;       the view object.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot::Draw, oImageRef=oImageRef, oViewRef=oViewRef, oWindowRef=oWindowRef

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  ;; Set up an error handler
  CATCH, theError
  IF theError NE 0 THEN BEGIN
    CATCH, /CANCEL
    ok = ERROR_MESSAGE(!Error_State.Msg + ' Returning...', TRACEBACK=1, /ERROR)
    RETURN
  ENDIF

  IF N_ELEMENTS(oImageRef) NE 0 THEN BEGIN
    ;; It's object graphics
    *oImageRef  -> SetProperty, DATA=*self.fractal
    *oWindowRef -> Draw, *oViewRef
  ENDIF ELSE BEGIN
    ;; It's direct graphics
    cgImage, *self.fractal, /Scale, /Keep_aspect
  ENDELSE

END


;+
; :Description:
;   The definition method for the DJP_Mandelbrot class.
;
; :Author:
;   David J Platten
;
; :History:
;   Change history::
;     Written, May 2013, DJP
;
; :Copyright:
;   Copyright (c) 2013, David J Platten
;-
PRO DJP_Mandelbrot__Define

  COMPILE_OPT IDL2 ;; Force all integer variables to be LONG and not SHORT. Force [] for arrays

  struct = { DJP_Mandelbrot,       $ ;; The object class name
             fractal   : PTR_NEW(), $ ;; The INTEGER array representing the fractal
             indices   : PTR_NEW(), $ ;; The LONG array representing the indices of fractal
             zed       : PTR_NEW(), $ ;; The COMPLEX array representing the iterated values
             coords    : PTR_NEW(), $ ;; The COMPLEX array representing the space being examined
             imageSize : [0, 0]   , $ ;; The real and imaginary pixel side-lengths of the arrays
             rMin      : 0.0D     , $ ;; The minimum real value
             rMax      : 0.0D     , $ ;; The maximum real value
             iMin      : 0.0D     , $ ;; The minimum imaginary value
             iMax      : 0.0D     , $ ;; The maximum imaginary value
             numIter   : 0        , $ ;; The number of iterations carried out
             extra     : PTR_NEW() }  ;; Extra keywords
END